Homework#2

Задание 1

lifexp <- readRDS("data/life_expectancy_data.RDS")

Задание 2

plot_ly(data = lifexp %>% 
          mutate(lbl = sprintf("Unemployment: %.2f Life expectancy: %.1f",
                               Unemployment, `Life expectancy`)),
        x = ~Unemployment, y = ~`Life expectancy`, 
        type = 'scatter', mode = "markers", 
        color = ~continent,
        text = ~lbl, hoverinfo = "text",
        marker = list(size = 10), alpha = 0.5) %>%
  layout(xaxis = list(title = "Unemployment, %"),
         yaxis = list(title = "Life expectancy, year", zeroline = FALSE),
         title = "Life expectancy vs. Unemployment, by Continent")

Задание 3

lifexp1 <- lifexp %>% 
  filter(continent %in% c("Africa", "Americas")) 

stat.test <- lifexp1 %>% 
  rstatix::t_test(`Life expectancy` ~ continent) %>% 
  add_xy_position(x = "continent")

ggboxplot(
  lifexp1, 
  x = "continent", y = "Life expectancy", 
  ylab = "Life expectancy, year", xlab = "Continent", 
  add = "jitter"
  ) + 
  labs(subtitle = get_test_label(stat.test, detailed = TRUE)) + 
  stat_pvalue_manual(stat.test, tip.length = 0) 

Задание 4

lifexp2 <- lifexp %>% 
  select(is.integer | is.numeric, -Year)

cor_plot <- lifexp2 %>%
  select(everything()) %>%
  psych::corr.test(adjust = "BH")  

corrplot(corr = cor_plot$r,
         p.mat = cor_plot$p,
         method = "color",
         order = "hclust")

cor_plot2 <- ggpairs(lifexp2,
        title = 'Correlations in lifeExp dataset',progress = F) +
    theme_minimal() +
    scale_fill_manual(values = c('#69b3a2')) +
    scale_colour_manual(values = c('#69b3a2'))

cor_plot2

Задание 5

lifexp2 <- lifexp %>%
  column_to_rownames('Country')%>% 
  select(is.integer | is.numeric, -Year)

lifexp2_scaled <- scale(lifexp2) %>% round(2)
lifexp2_dist<- dist(lifexp2_scaled, 
                        method = "euclidean")

as.matrix(lifexp2_dist)[1:6,1:6]

##                     Afghanistan  Albania  Algeria   Angola Antigua and Barbuda
## Afghanistan            0.000000 7.605169 6.330893 4.411814            6.647481
## Albania                7.605169 0.000000 2.620801 7.916066            3.357946
## Algeria                6.330893 2.620801 0.000000 6.318647            4.347229
## Angola                 4.411814 7.916066 6.318647 0.000000            8.092558
## Antigua and Barbuda    6.647481 3.357946 4.347229 8.092558            0.000000
## Argentina              7.922399 3.625245 3.459422 7.159518            4.967353
##                     Argentina
## Afghanistan          7.922399
## Albania              3.625245
## Algeria              3.459422
## Angola               7.159518
## Antigua and Barbuda  4.967353
## Argentina            0.000000

lifexp2_hc_ward <- hclust(d = lifexp2_dist, 
                        method = "ward.D2")

fviz_dend(lifexp2_hc_ward, 
          cex = 0.4, main = "Dendrogram, All variables numeric, Ward's method",
          k=5, k_colors = c("#cc0337", "#010d85", "#37953f", "#f98866", "#5BC7F2"),
          color_labels_by_k = TRUE,
          rect = TRUE)

fviz_dend(lifexp2_hc_ward, 
          cex = 0.4, main = "Dendrogram, All variables numeric, Ward's method(circular)",
          k=5, k_colors = c("#cc0337", "#010d85", "#37953f", "#f98866", "#5BC7F2"),
          color_labels_by_k = TRUE,
          type = "circular")

Задание 6

pheatmap(lifexp2_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = lifexp2_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5,
         cutree_cols = length(colnames(lifexp2_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

Строки объединены в 5 кластеров (исходя из дендрограммы выше). Довольно ярко выражена однородность между переменными GDP и GNI (для всех кластеров). Для первого и второго кластеров подобная однородность прослеживается для переменной Per Capita. Весьма примечателен второй клатер, включающий в себя узкую группу наблюдений. Также более менее однородным является первый и второй кластера для переменнной Hospital beds. Для переменных Measless Immunization, DPT Immunization и HepB3 Immunization ситуация аналогична предыдущей (особенно первый кластер).

Задание 7

lifexp_full.pca <- prcomp(lifexp2, 
                        scale = T,
                        ) 

summary(lifexp_full.pca)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 3.394e-16
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00

fviz_eig(lifexp_full.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(lifexp_full.pca, select.var = list(contrib = 5), col.var = "contrib")

fviz_contrib(lifexp_full.pca, choice = "var", axes = 1, top = 20) # 1

fviz_contrib(lifexp_full.pca, choice = "var", axes = 2, top = 20) # 2

fviz_contrib(lifexp_full.pca, choice = "var", axes = 3, top = 20) # 3

Наибольший вклад в изменении вариативности наблюдений вносит первая компонента PC1 (39.9%). Далее представлен PCA-plot, для которого число наиболее важных переменных для двух главных компонент сокращено до 5 и выделен вклад каждой из них (посредством длины стрелки переменной и её направления, всего 3 группы), при этом переменная Infant Mortality направлена обратно пропорциональна Immunization переменным. Также были отдельно отображены конкретные переменные, вклад которых наиболее существенен для первых 3 главных компонент (особенно приятен взору разбор третий компоненты, где первые две переменные объясняют львиную долю её вариации).

Задание 8

plott <- ggbiplot(lifexp_full.pca, 
         scale=0, 
         groups = as.factor(lifexp$continent),
         ellipse = T,
         alpha = 0.2,
         labels = rownames(lifexp2), 
         labels.size = 2.1) +
  theme_minimal()

plott

ggplotly(plott)

Задание 9

Можно резюмировать, что алгоритм PCA на этих данных сошёлся довольно скверно, т.к. две первые главные компоненты объясняют чуть более 50% изменений в данных - далее прирост становится и вовсе примерно одинаковым, постепенно уменьшаясь (смотри на строку Cumulative Proportion в таблице Importance of components) - можно предположить, что в данных как таковых корреляций между переменными и скрытых факторных переменных нет ==> потенциально доступен для использования регрессионный анализ.

Задание 10

umap_prep <- recipe(~., data = lifexp2) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors()) %>%  
  prep() %>%  
  juice() 

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(color = lifexp$continent), 
             alpha = 0.7, size = 2) +
  labs(color = NULL) 

Алгоритмы PCA и UMAP отработали схожим образом, дивергируя наблюдения по континентам приблизительно в одной манере - например, континенты Африка и Европа в обоих случаях находятся по разные стороны (в сравнении с графиком из задания 8) шкалы UMAP1 (аналог PCA1).

Задание 11

pca_run <- function(cols) {
  
  dat <- lifexp2 %>% select(-cols)

  lifexp_full.pca <- prcomp(dat, 
                          scale = T) 
  
  print(summary(lifexp_full.pca))
  
  print(fviz_eig(lifexp_full.pca, addlabels = T, ylim = c(0, 40)))
  
  print(fviz_pca_var(lifexp_full.pca, select.var = list(contrib = 5), col.var = "contrib"))
  
  print(fviz_contrib(lifexp_full.pca, choice = "var", axes = 1, top = 20)) # 1
  print(fviz_contrib(lifexp_full.pca, choice = "var", axes = 2, top = 20)) # 2
  print(fviz_contrib(lifexp_full.pca, choice = "var", axes = 3, top = 20)) # 3
  
  plott <- ggbiplot(lifexp_full.pca, 
         scale=0, 
         groups = as.factor(lifexp$continent),
         ellipse = T,
         alpha = 0.2) +
  theme_minimal()

  print(plott)

}

vec1 <- c("Life expectancy", "GDP", "Per Capita", "Hospital beds", "Sucide Rate")
pca_run(vec1) 

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4860 1.4445 1.06741 1.02232 0.90437 0.89062 0.79027
## Proportion of Variance 0.4415 0.1490 0.08138 0.07465 0.05842 0.05666 0.04461
## Cumulative Proportion  0.4415 0.5905 0.67189 0.74654 0.80496 0.86162 0.90623
##                            PC8     PC9    PC10    PC11    PC12   PC13      PC14
## Standard deviation     0.70284 0.57430 0.43406 0.36112 0.35414 0.2116 2.795e-16
## Proportion of Variance 0.03528 0.02356 0.01346 0.00931 0.00896 0.0032 0.000e+00
## Cumulative Proportion  0.94151 0.96507 0.97853 0.98785 0.99680 1.0000 1.000e+00

vec2 <- c("Unemployment", "Clean fuels and cooking technologies", "Measles Immunization", "Basic sanitation services", "Rural population")
pca_run(vec2) 

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.2979 1.3973 1.2549 1.13077 0.92041 0.87450 0.8111
## Proportion of Variance 0.3772 0.1394 0.1125 0.09133 0.06051 0.05463 0.0470
## Cumulative Proportion  0.3772 0.5166 0.6291 0.72042 0.78093 0.83556 0.8825
##                            PC8     PC9    PC10   PC11    PC12    PC13    PC14
## Standard deviation     0.69911 0.66898 0.55976 0.5061 0.30325 0.20399 0.06989
## Proportion of Variance 0.03491 0.03197 0.02238 0.0183 0.00657 0.00297 0.00035
## Cumulative Proportion  0.91746 0.94943 0.97181 0.9901 0.99668 0.99965 1.00000

vec3 <- c("GNI", "Mortality caused by road traffic injury", "DPT Immunization", "Tuberculosis treatment", "Non-communicable Mortality", "Rural population")
pca_run(vec3) 

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3544 1.2394 1.13818 1.01794 0.94990 0.87212 0.75704
## Proportion of Variance 0.4264 0.1182 0.09965 0.07971 0.06941 0.05851 0.04409
## Cumulative Proportion  0.4264 0.5446 0.64422 0.72392 0.79333 0.85184 0.89592
##                           PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.6480 0.60632 0.47451 0.38854 0.33849 0.27343
## Proportion of Variance 0.0323 0.02828 0.01732 0.01161 0.00881 0.00575
## Cumulative Proportion  0.9282 0.95650 0.97382 0.98544 0.99425 1.00000

В целом мы можем наблюдать скачки во вкладе тех или иных компонент в объяснении вариативности данных (что видно как по графику 1, так и по таблице Importance of components), но везде первая главная компонента остаётся самой значимой в процентном отношении. На PCI-плотах могут становиться значимые другие переменные после удаления случайных. Понятно, что также будет отличаться вклад тех или иных переменных в первый, второй и третий компоненты. Что касается biplot’ов тут можно лицезреть определённое единнобразие в расположение эллипсов, группирующих наблюдения по континентам (исключая набор переменных-стрелок, конечно).

Задание 12

lifexp3 <- lifexp %>% 
  mutate(african_flag = if_else(continent == "Africa", 1, 0),
         oceania_flag = if_else(continent == "Oceania", 1, 0)) %>% 
  select(is.integer | is.numeric, -Year)

lifexp_full.pca <- prcomp(lifexp3, 
                        scale = T) 

summary(lifexp_full.pca)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion  0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion  0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
##                           PC15    PC16    PC17    PC18    PC19    PC20
## Standard deviation     0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion  0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
##                             PC21
## Standard deviation     3.985e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

fviz_eig(lifexp_full.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(lifexp_full.pca, select.var = list(contrib = 5), col.var = "contrib")

fviz_contrib(lifexp_full.pca, choice = "var", axes = 1, top = 20) # 1

fviz_contrib(lifexp_full.pca, choice = "var", axes = 2, top = 20) # 2

fviz_contrib(lifexp_full.pca, choice = "var", axes = 3, top = 20) # 3

plott <- ggbiplot(lifexp_full.pca, 
         scale=0, 
         groups = as.factor(lifexp$continent),
         ellipse = T,
         alpha = 0.2) +
theme_minimal()

plott

Исходя из полученных результатов, можно предположить, что добавление dummy-переменных не оказало существенного влияния на объяснение вариативности данных (но заметен определённый вклад dummy-переменной в PC1 - той что маркирует записи из африканского континента). Потенциально, имея ввиду тот факт, что переменная continent не числовая и не участвует в анализе PCA, внедрение подобных dummy-деривативов могло бы улучшить анализ. Однако, с другой стороны, эти dummy переменные в какой-то степени дублируют уже имеющуюся переменную continent, которая используется для группировки наблюдений. При этом надо заметить, что такие dummy содержат меньше информации, чем оригинальные переменные.